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PARTICLE-BASED MULTISCALE MODELING OF 
CALCIUM PUFF DYNAMICS* 

ULRICH DOBRAMYSL 1 , STEN RUDIGER 2 , AND RADEK ERBAN 1 

Abstract. Intracellular calcium is regulated in part by the release of Ca 2+ ions from the endo¬ 
plasmic reticulum via inositol-4,5-triphosphate receptor (IP3R) channels (among other possibilities 
such as RyR and L-type calcium channels). The resulting dynamics are highly diverse, lead to local 
calcium “puffs” as well as global waves propagating through cells, as observed in Xenopus oocytes, 
neurons, and other cell types. Local fluctuations in the number of calcium ions play a crucial role 
in the onset of these features. Previous modeling studies of calcium puff dynamics stemming from 
IP3R channels have predominantly focused on stochastic channel models coupled to deterministic 
diffusion of ions, thereby neglecting local fluctuations of the ion number. Tracking of individual ions 
is computationally difficult due to the scale separation in the Ca 2+ concentration when channels are 
in the open or closed states. In this paper, a spatial multiscale model for investigating of the dynam¬ 
ics of puffs is presented. It couples Brownian motion (diffusion) of ions with a stochastic channel 
gating model. The model is used to analyze calcium puff statistics. Concentration time traces as 
well as channel state information are studied. We identify the regime in which puffs can be found 
and develop a mean-field theory to extract the boundary of this regime. Puffs are only possible when 
the time scale of channel inhibition is sufficiently large. Implications for the understanding of puff 
generation and termination are discussed. 
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1. Introduction. Intracellular calcium plays a major role in many signaling 
pathways, and regulates enzymatic activity, gene expression [7], and neural activity 
[60, 36]. In order to control a variety of cell functions, cells control the local cyto¬ 
plasmic calcium concentration via the exchange of Ca 2+ ions with the extracellular 
space and the release of Ca 2+ ions through channels situated on the membrane of 
reservoirs, such as the endoplasmic and the sarcoplasmic reticulum. Modulation of 
the Ca 2+ release regulates muscle contraction, pathway cross-talk and mitochondrial 
activity, and disruption of these processes is associated with various diseases such as 
early-onset Alzheimer’s [5, 57], heart failure [3, 2] as well as psychological conditions 
such as bipolar disorder and schizophrenia [6]. Hence, detailed knowledge about the 
underlying processes governing this signaling mechanism is required to allow progress 
in our understanding of these diseases. 

In this paper, we focus on the release of calcium ions from the endoplasmic retic¬ 
ulum (ER) via inositol-4,5-triphosphate receptor (IP3R) channels. Upon binding of 
Ca 2+ to binding sites on its cytosolic part, a channel opens and calcium ions flow from 
the ER into the cytoplasm. Channels occur in clusters of 10 to 20 channels [12, 47]. 
The opening of a single channel usually triggers the release of Ca 2+ from other chan¬ 
nels in the same cluster due to the increased Ca 2+ concentration in their vicinity [45]. 
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This mechanism results in a highly localized increase of the cytosolic calcium concen¬ 
tration. These “puffs” of calcium ions have been detected and analyzed in experiments 
by using fluorescent calcium buffers [58, 23]. 

Traditional modeling approaches of intracellular calcium dynamics are based on 
deterministic macroscopic rate equations [27, 11], however the intrinsically random, 
erratic nature of calcium signals in many cells or cell domains necessitates an ap¬ 
proach going beyond the deterministic regime [59, 24, 31, 44]. Progress has been 
made by recognizing the importance of number fluctuations in the binding to the 
channels [48, 18] and using hybrid models, where the deterministic calcium concen¬ 
tration is coupled to stochastic channel binding models [42]. Recently, it was found 
that local fluctuations stemming from diffusive noise (i.e. noise originating from the 
random movement of ions) have a crucial influence in calcium dynamics for clusters 
of intracellular channels, particularly in the inter-puff waiting time [ 22 ] but also on 
single-channel equilibrium behavior [56]. Stochastic effects were also investigated for 
L-type calcium channels and RyR channels in the dyadic cleft [34, 49, 28]. Other 
studies consider a Langevin equation governing the fraction of open channels in a 
cluster [54]. However, tracking the exact diffusion of individual ions in the complete 
computational domain is computationally intensive. 

Here, for the first time, we apply spatial stochastic multiscale methods to model 
the dynamics of calcium puffs including the release of Ca 2+ from IP 3 R channels and 
track individual ion positions in order to accurately incorporate diffusive noise. We 
take into account both activating and inhibitory channel properties. We study the 
dynamics of the calcium concentration as a function of the ion binding affinities and 
explore the regime where there exist puffs, as well as the regime in which channels 
do not close after their first opening. The paper is organized as follows: In the next 
Section 2 we discuss the spatial stochastic model for diffusion and channel gating used 
throughout this study. Section 3 discusses the multiscale approach we employ in order 
to reduce the computational effort required to track single ions and hence make this 
study feasible. In Section 4 we present results on the statistics of puffs extracted from 
simulated time series data, and study the transition between perpetually-open channel 
clusters and the parameter regime in which puffs can be observed. In addition, we 
develop a mean-field model for channel dynamics and use it to extract the boundary 
between the regimes. Finally, we summarize our findings in Section 5. 

2. Spatial stochastic model for intracellular Ca 2+ release. The spatial 
extent of our computational model consists of the three-dimensional domain £2 which 
models a part of the intracellular space. Ca 2+ ions are able to undergo free diffusion 
in Cl. They bind to and dissociate from binding sites on the channels, which are 
positioned on a small area of the domain boundary, corresponding to the membrane 
of the ER. The domain geometry and boundary conditions are specified in Section 2 . 2 . 
In the following sub-sections we discuss the components of this model in detail. The 
parameter values used throughout this study can be found in Table 1. 

2.1. Diffusion - Brownian dynamics. A versatile method for the simulation 
of particles in a solution is given by Brownian Dynamics (BD). Collisions of particles 
with solution molecules lead to overdamped dynamics and random forcing on a suffi¬ 
ciently long time scale [15]. Assuming that there are Q(t ) free ions in the simulation 
domain at time £, the equation for Brownian motion of ions is given by 

( 2 . 1 ) dX, = V2D dWj j = 1 , 2 ,..., Q(t) , 
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Parameter description 

Name 

Value 

Parameter values from the literature 

Free Ca 2+ diffusion constant [1] 

D 

220 pm 2 s -1 

Cytoplasmic Ca 2+ concentration [43] 

co 

0.02 pM 

Edge length of computational domain [43] 

L 

5 pm 

Rate of binding to activating site [43] 

Cta 

100 pM — 1 s -1 

Rate of unbinding from activating site [43] 

K 

20 s- 1 

Open channel current [8, 47, 51] 

Ic 

0.1 pA 

Number of channels in cluster [43, 12, 48] 

c 

9 

Chosen simulation parameters 

Spacing between channels in cluster 

t 

0.15 pm 

Rate of binding to inhibitory site 

Cti 

varies 

Rate of unbinding from inhibitory site 

hi 

varies 

Binding radius 

Q 

0.03 pm 

Unbinding radius 

a 

0.015 pm 

BD time step 

At 

0.1 ms 

Edge length of BD regime 

Lbd 

1 pm 

Compartment size 

h 

0.2 pm 


Table 1 

Parameter values for the spatial stochastic simulations, for physiologically relevant conditions [43]. 


where D is the diffusion coefficient of ions in solution, Xj (t) E D describes the tra¬ 
jectory of the j -th ion, and W j is a three-dimensional vector of independent Wiener 
processes. This approach dramatically reduces the dimensionality of the problem 
compared to molecular dynamics approaches wherein the degrees of freedom of every 
participating molecule need to be taken into account. Nevertheless, the computa¬ 
tional load is still high compared to deterministic PDE-based approaches to diffusion. 
There are a number of approaches for simulating (2.1) in the literature, ranging from 
discretization with a fixed time step [4] to event-based methods [52, 37]. In this paper, 
we discretize time with the time step At and use the Euler-Maruyama discretization 
of equation (2.1), i.e. the position of the j-th ion is updated according to 

(2.2) X j (t + At) = X j (t) + y/2DAt$ j , j m 1 , 2 ,..., Q(t) , 

where is a vector of three independent normally distributed random numbers with 
zero mean and unit variance. 

2.2. Domain geometry and boundary conditions. In our simulations, the 
computational domain is given as cube D = [0,L] 3 where the value of L is specified 
together with other parameters in Table 1. Ca 2+ ions are able to undergo free diffusion 
in D which we simulate using (2.2). A cluster of nine IP3R channels is positioned in 
a 3 x 3 grid with grid constant i = 0.15 pm, centered in the z = 0 plane, i.e. the 
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positions of nine channels in the cluster are given as 


(2.3) 
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We found no significant effects from varying the channel spacing £, hence we chose this 
particular cluster configuration for ease of implementation. Ions bind to and dissociate 
from binding sites on the channels. The boundaries of the computational domain at 
x = 0, x = L, y = 0, y = L, and z = L are const ant-concentration boundaries, 
hence they absorb and introduce ions, such that in the absence of open channels the 
concentration of Ca 2+ in the computational domain is held at its equilibrium value, 
Co, on average. The boundary at z = 0 is reflective, corresponding to the membrane 
of the ER. Hence, ions can enter and leave the domain via the boundaries, in addition 
to being introduced through open channels. 

If the channels are closed, then the average number of free ions in the computa¬ 
tional domain, (Q(£)), is equal to co\Q\ where \Q\ = L 3 is the volume of Q. Using 
our parameter values (see Table 1), cq\Q\ ~ 1.5 x 10 3 . The average number of ions 
in the simulation domain, (Q(t)), increases when the channels are open (by around 
two orders of magnitude). In order to simulate the system over sufficiently long time 
intervals, we will use a multiscale approach, which we describe in detail in Section 3. 

2.3. Stochastic channel binding model. The conformational changes be¬ 
tween the open and closed states of IP3R channels are controlled by the binding 
of Ca 2+ to activating and inhibitory binding sites. These channels consist of 4 sub¬ 
units, that can be itself in an active or a neutral/inhibited state. Each subunit has 
three different binding sites: An activating binding site for Ca 2+ ions, an inhibitory 
binding site for Ca 2+ ions, and an IP3 binding site. 

To accurately capture channel gating events, we employ a simplified DeYoung- 
Keizer model [11]. Here, we disregard IP 3 dynamics and consider the effects of IP 3 
only via their influence on the dissociation constant bi of Ca 2+ ions from inhibitory 
binding sites. Free Ca 2+ ions can bind to the activating and inhibitory sites, while 
bound ions can dissociate from occupied sites. Therefore, there are two reversible 
reactions in our model for each subunit of a channel: 


(2.4) 


Ca 2+ + {S a 


0} £ {S t 


b a 


1 } 


Ca 2+ + {S) = 0} 



{$ = !} 


The rates a a and eq describe the binding affinity of a Ca 2+ ion to an activating and 
inhibitory binding site, respectively, while the off-rates b a and bi describe the corre¬ 
sponding dissociation reaction rates. The variables S a and Si describe the binding 
site state and can take values of 0 and 1 . Channels consist of four subunits, with each 
subunit having one activating and one inhibitory Ca 2+ -binding site; see Figure 1 . A 
model subunit can then be in three distinct states: neutral (no binding site occupied), 
active (only the activating site occupied) and inhibited (if the inhibitory site is oc¬ 
cupied regardless of the state of the activating site). A channel then opens when at 
least three of its four subunits are in the active state. 

2.4. Ion binding dynamics. In order to precisely capture channel opening and 
closing dynamics, we need to implement a reversible BD binding model of Ca 2+ ions 
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Fig. 1. (a) Sketch of a channel, containing 4 subunits with an activating ion binding site (green 
square) and an inhibiting site (red circle), (b) State space of a single subunit. Ions can bind to the 
activating and the inhibiting sites. A subunit is active only if it has an occupied activating site as 
well as an empty inhibiting site. 


to their corresponding binding sites. Ions become binding candidates as soon as they 
enter a half-sphere of radius g around a binding site (the channels are positioned 
on the ER membrane, hence only the half-sphere above them is available for ion 
binding). They are then allowed to bind to the site with a probability P\ per time 
step that they spend in the binding region [35, 16]. An ion bound to an activating 
site (resp. inhibitory site) is allowed to dissociate with a probability 1 — exp(— b a At) 
(resp. 1 — exp (—bi At)) and placed at a distance of a (unbinding radius) from the 
binding site. We use the values of binding and unbinding radii, g and cr, as given in 
Table 1. Their values are chosen to be reasonable, such that no overlaps occur between 
channels and that they are large enough such that the BD simulation time step can 
be chosen reasonably large. Then we can pre-calculate the remaining parameter P\ 
(binding probability) before the start of simulations using the approach described 
in [35, Section 5]. 

2.5. Channel opening and calcium flux. Each channel has four subunits. 
Let S(f k (t) (resp. Sj :k (t))^ j = 1, 2 ,..., 9, k = 1, 2, 3,4, be the state of the activating 
(resp. inhibiting) site of the k -th subunit of the j -th channel in the cluster. Then the 
number of active subunits of the j -th channel at time t is 

k=1 

If (t) > 3, then the channel is considered to be in an open state. When a channel is 
open, new ions are introduced with a rate of 3.12 x 10 5 s -1 (corresponding to channel 
current Ic = 0.1 pA [8, 47, 51]) at the channel site, simulating the flux of Ca 2+ ions 
out of an active channel. The positions of released ions then evolve according to 
equation (2.2). 

3. Multiscale approach. During a puff event, a large number of ions are re¬ 
leased into the cytoplasm. If all channels are open continuously, then one can estimate, 
for the parameter values given in Table 1, that our computational domain may con¬ 
tain of the order of 10 5 ions during the peak of a puff. Since we are interested in 
time scales of minutes, tracking the individual position of this number of ions via 
Brownian dynamics becomes infeasible. However, high accuracy and individual ion 
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Fig. 2. Sketch of the computational domain It. It consists of two regions: The compartment- 
based part Vt \ Qbd (everywhere except for the red box in the bottom center, with the compartment 
size h illustrated in the left bottom of the picture) and the BD domain Qbd (shown as a red box in 
the bottom center). The boundary in both x- and y-directions as well as the boundary at z = L are 
constant-concentration boundaries, while the boundary at z = 0 is reflective. The channel cluster is 
positioned at the center of the z = 0 boundary (indicated as black squares), see (2.3). 


positions are only needed in the vicinity of channel sites in order to ensure accurate 
implementation of BD binding dynamics, described in Section 2.4. Therefore, we split 
our computational domain into two regions: a cube 


£Ibd — 


L — Lbd L + Lbd 


L — Lbd L + Lbd 


x [ 0 ,L B d\ 


containing the channel sites (red region in Figure 2), as well as the remaining space 
LL\Qbd in which we will use a coarser description of ion movement as described in 
the next subsection. Here, Lbd < L is the length of the edge of the cube LLbd- We 
use Lbd = L/b in our simulations (see Table 1). In particular, we use BD simulations 
in a small fraction of 1/5 3 ~ 0.8% of the computational domain Cl. 

3.1. Compartment-based model for diffusion. We subdivide the region Cl \ 
LLbd into compartments (cubes) with size h (illustrated in the bottom left of Figure 2) 
and employ a compartment-based method to simulate the movement of ions [17]. 
Ions are allowed to move between adjacent compartments with a rate of d = D/ti 2 . 
The compartment-based algorithm only stores and evolves the number of ions in 
each compartment (rather than following individual ions). An event-based stochastic 
simulation algorithm is used to efficiently simulate this system. Several equivalent 
methods have been developed in the literature, such as the Gillespie algorithm [26], 
the Next Reaction Method [25], the Next Subvolume Method [14, 30], as well as the 
Optimized Direct Method [9]. 

In this paper, we employ the Next Reaction Method [25]. For each possible move 
between two neighboring compartments, a putative time for the next jump of an ion 
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to occur is calculated. It is given by 


(3-1) 


ln(r) 


dA c (t ) 


where t is the current simulation time, d = D/h 2 is the jump rate per one ion, A c (t ) is 
the current number of ions in the compartment from which an ion is jumping and r is a 
random number uniformly distributed in the interval (0,1). Clearly, the putative jump 
time (3.1) is infinity if A c (t) = 0, i.e. if the corresponding compartment is empty. The 
putative times are smaller (on average) if the corresponding compartment contains 
more ions. 

The putative jump times are inserted into a priority queue (a heap data struc¬ 
ture [10]), which enables us to efficiently extract the earliest jump time and thus 
the next jump. This move is then performed, and the numbers of ions in compart¬ 
ments (and the corresponding putative times and their entries in the priority queue) 
are updated. We then iterate this process by finding the minimal putative time and 
performing the corresponding ion jump at each iteration. 

At the boundaries of the computational domain D (except for the boundary at 
z = 0, which is reflective), ions are absorbed (i.e. they leave the domain) or can 
enter with rates consistent with a constant equilibrium concentration of Co outside 
the domain. The jump rate from outside the domain into a compartment just inside 
the domain boundary is cqDK which is used, instead of dAc{t ), in (3.1) to compute 
the corresponding putative times. 

3.2. Coupling the BD simulation in Qbd with the compartment-based 
approach in Q\Qbd» Several methods exist in order to couple the BD and compart¬ 
ment-based methods across their interface, such as the two-regime method (TRM) [19, 
20] or the ghost-cell method [21]. Here, we employ the TRM. In order to accurately 
capture diffusion across the interface, the jump rates from adjacent compartments 
into VLbd need to be adjusted from the bulk rate d = D/h 2 to the interface rate [19] 



(3.2) 


These jump rates are used, instead of d, in (3.1) to compute the corresponding putative 
times. If the chosen move in the compartment regime is a jump from a compartment 
adjacent to the interface into the occupancy number of the compartment is 

reduced by one and a new ion is added in Qbd at the distance x from the interface 
which is sampled from the probability distribution [19] 



(3.3) 


where erfc is the complementary error function. 

The rate (3.2) and the distribution (3.3) are used to transfer ions from D \ VLbd 
into £Ibd- In the opposite direction, the TRM transfers any ion from Qbd which 
during the time step interacts with the interface. For a detailed discussion and the 
derivation of the probabilities and jump rates above, please see references [19, 20, 21]. 
We employ the particle-based simulation library package Tyche which implements the 
TRM [38]. The TRM has also been recently implemented in Smoldyn [39]. 
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Fig. 3. Example concentration time trace for a^ = 1 pM 1 s 1 and bi = Is 1 . 

4. Results. In the following section we present the results of our findings. Sim¬ 
ulations were performed with the parameter values listed in Table 1 unless noted 
differently. In Section 4.1, we discuss Ca 2+ puff statistics from simulation runs. We 
follow this by investigating the regimes in which puffs are visible in Section 4.2. A 
simple mean-field theory is then presented in Section 4.3. It describes the transition 
between puffs and perpetually open channels. 

4.1. Puff statistics. Figure 3 displays a sample time trace of the Ca 2+ con¬ 
centration in the computational domain for a* = lpM -1 s -1 and bi = Is -1 for a 
time range of 20 s. The ion concentrations shows erratic calcium puffs with an am¬ 
plitude of up to 0.5 pM. The constant background of Co = 0.02 pM corresponds to 
approximately 12 ions per pm 1 * 3 * S and justifies our use of Brownian dynamics because 
fluctuations in such a small number of ions influence binding to the channels. The 
puffs are characterized by a sharp, almost instant increase in concentration, followed 
by an approximately exponential decay after the channel cluster closes. Puffs are sep¬ 
arated by a refractory phase in which channel subunits are inhibited and the cluster 
cannot open. 

We analyze the concentration data given as the time values tj and concentration 
values Cj , j = 1,2,..., TV, (where TV is the index of the last data point of the simulation 
run) by choosing a puff-starting threshold concentration value 

(4.1) c t = (cj) + yVar[ C j], 

i.e. one standard deviation above the mean concentration (cj). For the data set shown 
below, the threshold was c t « 0.1 [iM. The averaging in (4.1) is taken over all values 
of j, j = 1, 2 ,..., TV, i.e. the sample mean and variance in (4.1) are estimated by 

1 N x N 

(4.2) ( Cj ) = — ^ Cj and Var[ Cj ] = ( Cj - ( Cj )). 

j=i j=i 

Puffs are identified by the concentration crossing the threshold point (4.1). Puff 
ending points are identified via a re-crossing of this value. The crossing points are 
given in the ordered index set 

S = {j | cj -1 < c t and Cj >c u j = 2 ,3,..., TV - 1} 

(puff starting indices) and 

£ = {j I ^ > c t and c j+1 < c t , j = 2,3,..., TV - 1} 
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(puff ending indices). The sets S and £ are enumerated by the puff index 1 < k < N p , 
where N p = \£\ is the number of finished puffs in the data set, i.e. we disregard the 
last index in S if the last puff did not finish and denote 

S = {si, s 2j • • • 5 s n p } and £ = {ei, e 2 ,..., e n p }, 

where the elements in S (resp. £) are ordered, i.e. s\ < s 2 < ••• < sn p (resp. 
e\ < e 2 < • • • < ejv p ). The inter-puff times are then given by 

(4.3) T = ji afc+1 - t ek | s k+ 1 e S, e k e £, 1 < k < N p - l| H (t p , oo), 

where the threshold t p — 0.25 s. This threshold filtering decreases the number of puffs 
considered by removing the puffs which are not well separated. Indeed, the inset in 
Figure 4(a) shows that the distribution of interpuff times below the threshold follows 
an exponential distribution and is due to random channel reopenings stemming from 
residual calcium. The distribution of inter-puff times (histogram of set T) is shown in 
Figure 4(a) for a set of 954 puff intervals. Its shape is similar to a Gamma distribution 


(4.4) 


p(r) 


_ T ^- 1 p- T / 6) 

r(«)0* 


Indeed, when choosing k = 2.82 and 0 = 0.64 s in equation (4.4) (such that the 
resulting distribution has the same mean and variance as our data), the comparison 
is good. As the Erlang distribution (the special case of a Gamma distribution with 
k G N) results from summing k exponentially-distributed random variables, we can 
speculate that an effective model of calcium channel opening might be described 
by a three-stage binding/unbinding process. The maximum of the inter-puff time 
distribution, r m « 1.81 s, can be interpreted as the typical time between puffs in this 
system. 

The distribution of puff amplitudes are determined by taking the maximum value 
in a previously-determined puff interval 

A= \ a k \a k = max Cj, where Sk G 5, G £, 1 < k < N p 


The amplitude distribution (histogram of set A) is shown in Figure 4(b) and mirrors 
the broad distribution characterized in experiments [13]. 

The puff durations are calculated by considering the full duration at half max¬ 
imum, which is a common criterion for determining the puff duration in the liter¬ 
ature [13]. To this end, we calculate the indices of crossings of the half maximum 
threshold on the rising slope of a puff 

ft = jnfe | r k = max{j | cj-i < a k /2 and cj > a k /2 and s k < j < e k }, 
where s k G <S, e k G £ , a k G A, 1 < k < N p j , 
and the falling slope of a puff 

J 7 = [fk | fk = min{j | Cj > a k /2 and c j+1 < a k /2 and s k < j < e k }, 
where s k G 5 , e k G £ , a k G A, 1 < k < N p > . 
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Interpuff time [s] 


(b) 



(c) 



Duration [s] 


(d) 



Amplitude A [/iM] 

Fig. 4. Calcium puff statistics for a set of N p = 1151 puffs from simulations with a total 
duration of 3 simulation hours. The mean and standard deviations of the distributions are indicated 
by the solid and dashed red lines, respectively, (a) Inter-puff time distribution for 545 puff intervals 
after threshold filtering (4.3). The mean inter-puff time is 1.81 d= 1.077 s. The green dashed line 
displays Gamma distribution (4.4) with n = 2.82 and 0 = 0.64 s. The inset shows the distribution 
of interpuff times smaller than 250 ms with a mean of 0.59 d= 0.061s. (b) Amplitude distribution. 
The mean puff amplitude is 0.34 =b 0.12pM (indicated by the red vertical lines), (c) Puff duration 
(full duration at half amplitude) distribution. The mean puff duration is 0.04 =b 0.02 s. (d) A two- 
dimensional histogram showing the strong correlation between the average number of open channels 
(over a puff’s duration) with the puff amplitude. The black line indicates a linear fit to the correlation 
data. 
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The puff durations are then calculated via 


= {tf k - t r J f k G F, r k e TZ, l<k< N p }. 

The distribution of puff durations (histogram of set V) is plotted in Figure 4(c). 
The distribution is sharply peaked around 0.03 s and has a long tail stemming from 
repeated cluster re-openings during the concentration decay phase. The mean puff 
lifetime is within the range 40-70 ms measured in experiments [13]. 

Figure 4(d) displays the correlation between the number of open channels aver¬ 
aged over a puff’s duration M and the puff amplitude A. By finding the average 
M(A ) for each value of A, and performing a linear fit, we find the linear relationship 
M(A) = 2.75A + 0.23. 

4.2. Inhibition transition. Figure 5 shows the results of a simulation run with 
two different values of inhibitory site dissociation rate bi, illustrating two different 
modes of behavior. We use the same value of cq in all four panels (the rest of parameter 
values are given in Table 1). Figure 5(a-b) show the results of a simulation run 
with bi = 5 s -1 (corresponding to a dissociation constant of K& = bi/ai = 50 pM). 
These simulation parameters are consistent with the modeling of calcium puffs in 
the literature and used in hybrid PDE-based models [43]. Figure 5(a) shows the 
Ca 2+ concentration in the computational domain £2, while Figure 5(b) displays the 
fraction of occupied inhibitory sites in the channel cluster as a function of time. Due 
to the large dissociation rate, the necessary level of bound inhibitory sites can never 
be sustained long enough for all the channels in the cluster to close at the same 
time. Therefore, the Ca 2+ concentration is kept at a level where empty activating 
sites are immediately filled and the channel cluster stays perpetually open. Hence, 
puff termination in this case can be speculated to be facilitated by a mechanism 
other than channel inhibition, such as ER Ca 2+ reservoir depletion or a mechanism 
involving dissociation of IP 3 [40]. We checked the influence of lowering the binding 
radius g to 6 nm and found that it had no effect on the channel cluster closing. 

In contrast, a lower inhibitory site dissociation constant yields well-defined cal¬ 
cium puffs in our particle-based simulation scheme. Figure 5(c-d) display data from 
simulation runs with bi = 0.1s _1 (which corresponds to Kp = bi/di = lpM). Here, 
inhibitory sites binding is sustained on a high level for a long enough time such that all 
channels close and the excess Ca 2+ is removed. Hence there are time intervals when 
the cluster concentration decreases and reaches the equilibrium Ca 2+ concentration 
Co. Puffs are therefore clearly delineated and separated with a well-defined inter-puff 
time. 

Hence there exist two regimes: a puff regime [Figures 5(c) and 5(d)] and a regime 
with perpetually open channel clusters [Figures 5(a) and 5(b)]. We now proceed to 
characterize concentration and open channel time traces to find under which condi¬ 
tions well-defined puffs are possible. To this end, we use the “puff score” character¬ 
ization function introduced in [29], which quantifies the spike-ness of a given time 
trace of the number of open channels. We denote the number of open channels at 
time point tj by Oj, j = 1 , 2 ,..., TV, i.e. Oj E { 0 , 1 , 2 ,..., C} where C = 9 in our 
simulations. Then the puff score is defined by 

(4-5) [PS](a i ,6 i ) = ^ Y ^. 

where the averages are again take over all values of j, j = 1, 2,..., TV (compare with 


ll 





Time [s] 



Time [s] 


Fig. 5. (a) Calcium concentration and (b) fraction of inhibited channel subunits as a function 
of time for a\ = O.lpM -1 s —1 and bi = 5s —1 . Other parameter values are given in Table 1. Note 
that the channel cluster stays perpetually active in this configuration, (c) CVi 2+ concentration and 
(d) fraction of occupied inhibitory binding sites in the channel cluster as a function of time for 
the inhibitory site binding parameters ai = 0.1 pM -1 s -1 and bi = O.ls -1 . Other parameter values 
are given in Table 1. Note that the time scale of inhibition decay is long enough for the calcium 
concentration to decay completely, thus the channel cluster closes and we can identify well-defined 
puffs. 


12 





















|-| 0.50 

- 0.45 

- 0.40 

- 0.35 

- 0.30 ~ 

| rO 

- 0.25 J 

- 0.20 £ 

- 0.15 

- 0.10 
I 0.05 
™ 0.00 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
Inhibitory site binding a i [piM~ l s~ l ] 



Fig. 6. Phase plot of ai and bi. Color indicates the puff score (4.5) for the given set of 
parameters listed in Table 1. The black line shows the numerically-determined phase boundary from 
the mean-field model (4.6)-(4.8). 
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Fig. 7. Phase plot of ai and bi for the hybrid scheme described in the text. Color indicates the 
puff score (4.5) for the given set of parameters listed in Table 1. For comparison, the black line shows 
the numerically-determined phase boundary from the mean-field model (4.6)-(4.8) for Figure 6. 


(4.2)). The puff score (4.5) can take values in [0,1]. A puff score greater than 0.25 
indicates channel excitability and therefore the existence of puffs in the system. 

In order to visualize the two parameter regimes, we performed simulations for 
different inhibitory site binding parameters cq G [0.1,1] pM -1 s -1 and bi G [1,10] s -1 
for a simulation time duration of 100 s. We extracted the number of open channels 
over time and calculated the puff score (4.5). Figure 6 displays this quantity. The 
color indicates the value of [PS](a^,6i). The transition between the two regimes is 
not sharp, but gradual, especially for higher a^. This is due to prolonged channel re- 
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Fig. 8. Phase plot of a^ and bi with a low diffusion constant D = 20pm 2 s — 1 . The other 
parameters are given in Table 1. Color indicates the puff score (4.5). The black line shows the 
numerically-determined phase boundary from the mean-field model (4.6)-(4.8). Note that the scale 
of the y-axis here is different from Figures 6 and 7. 


openings becoming more likely due to faster dissociation of bound inhibiting ions when 
bi approaches the transition. The phase boundary is consistent with a dissociation 
constant of 4pM < Kp < 10 pM. 

To directly compare our results with previously-used schemes from the literature, 
we devised a simplified hybrid scheme: The main difference between our method and 
hybrid simulation algorithms is the much lower calcium concentration in the vicinity of 
open channels and the resulting weaker channel inhibition. Therefore, in the simplified 
hybrid scheme, whenever a channel is open, we do not use the particle-based binding 
described in section 2.4. Instead we assume a constant high calcium concentration 
of cb = 150 pM (consistent with results from hybrid simulations [43]) to generate 
random binding events to inhibitory sites with a rate of a^ce, irrespective of any ions 
in the vicinity. In all other respects, the simulation proceeds as described previously. 
Figure 7 shows the resulting map of puff scores. 

In order to study the influence of buffers (whose main effect is to slow down 
ion diffusion [33]) on the boundary of the puff regime, we performed a similar set 
of simulations with a lower Ca 2+ diffusion constant of D = 20 pm 2 s -1 . Because 
ions bound to buffer molecules can be viewed as slowly diffusing ions, lowering the 
diffusion constant is a way to model the effects of buffers [53] . The result is shown in 
Figure 8. With a lower diffusion constant, the calcium concentration in the channel 
cluster nanodomains decays more slowly. For puffs to exist, the time scale of the 
decay of inhibitory site binding needs to be longer than the time scale of Ca 2+ decay. 
Hence, the boundary separating the two regimes is pushed to smaller values of the 
inhibitory dissociation rate bi, corresponding to approximately 1 pM < Kb < 2pM. 
Hence, the effective diffusion constant plays an important role in determining the 
boundary between the two regimes. Note, that the overall variation of the puff score 
[PS](a$, bi) is smaller compared to the case of D = 220 pm 2 s -1 , therefore Ca 2+ puffs 
become less pronounced with slower ion diffusion. Ca 2+ buffers lead to additional 
extrinsic noise due to binding and unbinding of calcium ions to buffer molecules that 
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enhances Ca 2+ fluctuations which might have effects on the puff statistics on top of 
what we presented here [55]. 

4.3. Mean-field model. In order to find an approximate phase diagram to 
determine the parameter regimes in which calcium puffs occur, we develop a simplified 
non-spatial model. To this end, we consider an ensemble of identical channels that 
interact by a shared calcium domain (all-to-all coupling). The concentration c(t) in the 
1 pm 3 cube around a channel undergoes exponential decay with a phenomenological 
decay parameter A due to diffusive equilibration, and a linear increase with an ion 
influx rate v when the channel is open. The open states of the ensemble of channels is 
determined by the mean number of occupied activating and inhibitory binding sites 
per channel, a(t) and 6(t), in a similar way as in the spatial model above. In a simple 
representation of the subunit dynamics and their binding cooperativity we require 
that a channel is open at a given time t if a(t) >3 and b(t) < 2. Here the variables 
a(t) and b(t) describe how many subunits, on average, have activating and inhibitory 
ions bound to their respective binding sites. They evolve according to the mass-action 
rate equations corresponding to the reactions (2.4). Hence the model equations are: 

dc 

(4.6) — = 0(a — 3) 0(2 — b) v — A (c — Co) , 

(4.7) ^ =a a c(4-a) -b a a , 

(4.8) ^ = aic(4-b) - bib . 

Here, 0(x) is the Heaviside function with the properties 


©(a;) = 


0 

< 1/2 
1 


forx < 0 
for x = 0 
iorx > 0 


The first term on the right hand side of equation (4.6) describes the above-mentioned 
channel openings: The channels only open if three subunits are active and not in¬ 
hibited. The influx rate is determined via the channel current v = ( 2eV)~ 1 Ic = 
518.28 pMs -1 (where e = 1.602 x 10 -19 C is the electron charge; the in-flowing ions 
are assumed to be spread over a volume of V = 1 pm 3 ). This value is also consistent 
with influx rates extracted from the rising flanks of puffs in our simulations. The 
exponential decay parameter A was determined by fitting an exponential decay to cal¬ 
cium puff simulation data. The parameters of the mean-field model are summarized 
in Table 2. 

Figure 9 shows data from three representative numerical solutions of the ODE 
system (4.6)-(4.8) where the inhibitory binding rate is set to a* = 0.5 pM -1 s _1 . The 
initial conditions are c(0) = Co, a(0) = 4 and b( 0) = 0. The concentration c(t) in 
Figure 9(a) shows a puff with its characteristic exponential decay to cq. For bi = 2s -1 , 
c(t) returns to its equilibrium value Co, while for higher values, it oscillates at a 
high concentration level. These two states correspond to the “puff” regime and the 
“always open” regime, respectively. Figure 9(b) shows the corresponding trajectories 
in a(t)-b(t ) space. When the trajectory, after its initial transient phase, does not 
touch the line segment {(a, 2) for a E (3,4)}, highlighted in gray in Figure 9(b), the 
concentration will decay to Co and the system returns to its original state. If, however 
the a(t)-b(t ) trajectory hits this boundary, the system stays excited and the channels 
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Co 

0.02 pM 

Background Ca 2+ 

V 

5.18 x 10 2 pM s —1 

Source rate 

A 

22.9 s" 1 (D = 220 pm 2 s" 1 ) 
2.2 s" 1 (D = 20 pm 2 s" 1 ) 

Ca 2+ decay 

a a 

100 pM -1 s -1 

Rate of activating site binding 

K 

20 s" 1 

Rate of activating site unbinding 

CLi 

[0.1, 1] pM -1 S -1 

Rate of inhibiting site binding 

hi 

[0.1,8] s" 1 

Rate of inhibiting site unbinding 


Table 2 

Parameter values for the ODE mean-field model (4.6)-(4.8). 



Time t [s] 



Fig. 9. (a) Concentration over time and (b) phase space trajectories of the mean field system 

(4.6)-(4.8). The ratchet-like oscillation visible in the concentration traces for bi =4s — 1 and 5s —1 
stems from the channels being re-activated repeatedly due to insufficient inhibition. 


stay perpetually open. This can be translated into a temporal criterion by viewing 
the time evolution of the system as a two-step process: (1) Find the time r such 
that b(r) = 2 (i.e. the time until the second Heaviside function in equation (4.6) 
becomes zero); (2) find T\> t such that b{r\) = 2 and > r such that afo) = 3. 
If t 2 < Ti, well-defined puffs are possible, otherwise c(t) stays elevated and channels 
stay perpetually open. 

We now proceed to find the mean-field phase boundary between the “puff” regime 
and the “perpetually-open” regime. Steady-state analysis of equations (4.6)-(4.8) does 
not yield the correct results in the latter regime due to the discontinuous nature of 
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the Heaviside functions in equation (4.6) and the resulting temporal criterion above. 
The system rapidly switches between the inhibited state with b(t) > 2 and the non- 
inhibited state with b(t) < 2, leading to the observed decaying oscillations of c(t ) in 
Figure 9(a). Hence, we numerically integrate equations (4.6)-(4.8) to test if a given 
parameter set a^, bi falls into one of the two regimes. We let the system first evolve 
until b > 2 and continue until either a(t ) becomes smaller than 3 or b{t) becomes 
smaller than 2. We then identify the regime the system is in according to the criterion 
described above. A bisection algorithm is used to find the boundary bi{ai) between 
the two regimes. The black lines in Figures 6, 7 and 8 show the extracted phase 
boundaries for the parameter values given in Table 2. 

5. Conclusions. In this paper, we have reported a novel application of a particle- 
based spatial algorithm for diffusion to investigate the influence of diffusive noise on 
the dynamics of intracellular calcium release. Particle number noise in calcium micro¬ 
domains has attracted interest in recent studies [22, 55, 56, 32], also for L-type and 
RyR calcium channels [34, 49, 28], and it is important to clarify whether calcium 
diffusion as an additional noise source needs to be incorporated to obtain a better 
understanding of sub-cellular calcium signals. 

In order to make this study feasible, we split the domain into a compartment- 
based regime and a Brownian dynamics regime, coupled via the TRM. This allowed 
us to model the dynamics of full puffs, including release of a realistic number of ions 
and inhibition dynamics. We extracted concentration time traces and analyzed the 
resulting puff statistics. The inter-puff time distribution as well as the distributions 
of puff amplitudes and lifetimes agree qualitatively with experimental data in the 
literature [23]. 

We then proceeded to analyze the binding parameter regimes under which well- 
defined Ca 2+ puffs are possible. We found that, surprisingly, an inhibitory binding 
site dissociation constant Kd = 50 pM consistent with the literature and patch-clamp 
experiments [46], does not yield puffs in our model. In this parameter regime, channels 
stay perpetually open. In order to investigate the transition between well-defined 
puffs and perpetually-open channels, we characterized calcium concentration traces for 
various combinations of the inhibitory site binding parameters. The phase boundary 
visible in our data is consistent with a dissociation constant in the region 4 pM < Kp < 
10 pM. Lower values of the Ca 2+ diffusion constant yield a phase boundary at smaller 
inhibitory site dissociation rates and thus an even smaller dissociation constant. Given 
the reliable puff generation and termination in previous studies based on fitted gating 
models with large Kjj [43, 50] this is an unexpected result. 

Why does our model not show the same robust termination at large Kd as the 
hybrid approaches [41, 43, 44]? In the latter models the binding and unbinding to the 
receptors is stochastic but the calcium distribution is calculated from deterministic 
reaction-diffusion equations. To test for the differences of puffs in both models we 
have performed a set of simulations using large local Ca 2+ concentrations for open 
channels, similar to what is obtained in the hybrid method. These simulations confirm 
the robust termination in the hybrid scheme even at large Kp. In a hybrid approach, 
a deep inhibitory state is achieved owing to the large local nano-domain around each 
open channel displayed in the solution of deterministic calcium equations [41]. Thus, 
in the hybrid model, a channel is not inhibited by shared calcium in the domain but 
by its own released calcium. This effect has been termed self-inhibition [50]. 

In our simulations at large Kjj , however, there is insufficient inhibition to the 
channels during the early phase of a puff. This means that no or perhaps only one 
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subunit per channel binds inhibitory calcium, while reliable inhibition requires binding 
of three or four calcium ions. Our results suggest that diffusive noise mixes calcium in 
the cluster domain and diminishes the localized domains around open channels and 
self-inhibition. Thus we are led to a model of Ca 2+ puffs that is very different from the 
previous hybrid model. Our diffusive model allows inhibition only from a much less 
localized concentration profile and can therefore only be achieved at a much smaller 
dissociation constant of inhibitory binding sites. This conclusion is also supported by 
the the mean-field ODE model that we have devised and that captures the average 
binding state of the cluster’s activating and inhibitory sites as well as the resulting 
Ca 2+ concentration in a shared and well-mixed micro-domain. This model displays a 
sharp phase boundary between the two regimes, which agrees well with the data from 
our spatial simulations. 

Evidence for an inhomogeneous or homogeneous calcium distribution in the clus¬ 
ter is hard to obtain directly from our simulations because of the short lifetime of nano 
domains. In any case, our study highlights the role of the local calcium concentration 
in the termination of puffs and shows that puffs are very sensitive to fluctuations of 
residual calcium remaining after channel closing. It has to be noted though, that, 
apart from the diffusive noise, there are other differences of the current BD setup and 
the former hybrid approaches. Notably, these include the presence of calcium-binding 
buffers and the action of SERCA pump terms on the ER membrane boundary [43] and 
it remains to be analyzed to which extent these differences affect puff termination. 
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